Abstract: The goal of this study is to look at territory size of the Louisiana Waterthrush during breeding season at three locations; Wall Branch stream located in Rotary Park in Clarksville Tennessee 37040, Dry Creek located in Cheatham Wildlife Management Area in Ashland City Tennessee 37015, Big Hollow located in Beaman Park in Ashland City Tennessee 37015. Mapping these territories in these streams will show not only the territories of the Louisiana Waterthrush present but also possible areas of concentrated use that may aid in the location of their nesting sites.
We will be making KDE maps so we will need to register with google maps to use their satelite version of our map area. The link below is where a lot of the following code came from. https://cfss.uchicago.edu/notes/raster-maps-with-ggmap/ You must register for google maps to get an AI (I used the one year free trial for this and they don’t auto charge once the year is over). OpenMapSource is another option as well if your java supports this. Once you have registered create your base map using google maps address or coordinates.
register_google("AIzaSyDKCi-txHXwSLjB-UURgX1eLzhNVcsXMRg")
basemap<- get_map(location = "1539 Dry Creek Rd
Ashland City, TN 37015", zoom=15)
## Source : https://maps.googleapis.com/maps/api/staticmap?center=1539%20Dry%20Creek%20Rd%0AAshland%20City,%20TN%2037015&zoom=15&size=640x640&scale=2&maptype=terrain&language=en-EN&key=xxx-txHXwSLjB-UURgX1eLzhNVcsXMRg
## Source : https://maps.googleapis.com/maps/api/geocode/json?address=1539+Dry+Creek+Rd%0AAshland+City,+TN+37015&key=xxx-txHXwSLjB-UURgX1eLzhNVcsXMRg
ggmap(basemap)
Next, we can put our LOWA encounter cooridinates on the basemap we just created.
DCmap<-ggmap(basemap)+ geom_point(data=DryCreek, aes(Longitude, Latitude)) +
geom_point(data = DryCreek, aes(Longitude, Latitude), size = 0.1) +
theme(axis.title = element_text(face="bold")) + labs(x="Longitude", y="Latitude") +
theme(panel.border = element_rect(colour = "black", fill=NA, size=1)) + theme_bw()
plot(DCmap)
## Warning: Removed 1 rows containing missing values (geom_point).
## Warning: Removed 1 rows containing missing values (geom_point).
Now it’s time to create KDE and MCP, but first we have to format a few things. It will be important to first put your data into SpatialPoints we will be using CRS (coordinate reference system) which is the format for most GPS data. The link below is where a lot of the following code came from. https://mhallwor.github.io/_pages/activities_GenerateTerritories
DCpts <- sp::SpatialPoints(coords = cbind(DryCreek$Longitude, DryCreek$Latitude),
proj4string = sp::CRS("+init=epsg:4326"))
head(DCpts)
## class : SpatialPoints
## features : 1
## extent : -87.18876, -87.18876, 36.18876, 36.18876 (xmin, xmax, ymin, ymax)
## coord. ref. : +init=epsg:4326 +proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0
## put all the data into a single SpatialPointsDataFrame (spdf)
DC_spdf <- sp::SpatialPointsDataFrame(DCpts, DryCreek)
head(DC_spdf)
## # A tibble: 6 x 4
## Name Longitude Latitude Number
## <chr> <dbl> <dbl> <dbl>
## 1 Dry Creek -87.2 36.2 1
## 2 Dry Creek -87.1 36.2 2
## 3 Dry Creek -87.1 36.2 3
## 4 Dry Creek -87.1 36.2 4
## 5 Dry Creek -87.1 36.2 5
## 6 Dry Creek -87.1 36.2 6
#although we don't have seperate species we need to let R know that all points are for one stream system and this will create a list of one
DC_sep <- split(x = DC_spdf, f = DC_spdf$Name, drop = FALSE)
Time to make the MCP.
DCmcp <- lapply(DC_sep, FUN = function(x){rgeos::gConvexHull(x)})
##this makes polygon from the list of one when we told R that all variables are from one stream system
DCmcp <- mapply(DCmcp, names(DCmcp),
SIMPLIFY = FALSE,
FUN = function(x,y){x@polygons[[1]]@ID <- y
return(x)})
DCmcp <- do.call(rbind,DCmcp)
DCmcp <- SpatialPolygonsDataFrame(Sr = DCmcp,
data = data.frame(Bird = names(DCmcp)),
match.ID = FALSE)
plot(DCmcp)
Now that we have an MCP for the stream, onto KDE.
## Step one: do least squares cross-validation to estimate bandwidth (you may get a warning message but keep going)
bw <- lapply(DC_sep, FUN = function(x){ks::Hlscv(x@coords)})
## Warning in ks::Hlscv(x@coords): Data contain duplicated values: LSCV is not
## well-behaved in this case
## Step two: generate kde
DC_kde <-mapply(DC_sep,bw,
SIMPLIFY = FALSE,
FUN = function(x,y){
raster(kde(x@coords,h=y))})
# This code makes a custom function called getContour.
# Inputs:
# kde = kernel density estimate
# prob = probabily - default is 0.95
getContour <- function(kde, prob = 0.95){
# set all values 0 to NA
kde[kde == 0]<-NA
# create a vector of raster values
kde_values <- raster::getValues(kde)
# sort values
sortedValues <- sort(kde_values[!is.na(kde_values)],decreasing = TRUE)
# find cumulative sum up to ith location
sums <- cumsum(as.numeric(sortedValues))
# binary response is value in the probabily zone or not
p <- sum(sums <= prob * sums[length(sums)])
# Set values in raster to 1 or 0
kdeprob <- raster::setValues(kde, kde_values >= sortedValues[p])
# return new kde
return(kdeprob)}
DC_95kde <- lapply(DC_kde,
FUN = getContour,prob = 0.95)
These next plots put MCP on top of KDE
plot(DC_kde[[1]])+
plot(DCmcp[1,],add = TRUE)
## integer(0)
Time to overlap KDE onto the base map. Code found from http://data-analytics.net/cep/Schedule_files/geospatial.html
KDEmap<-DCmap +
stat_density2d(aes(x = Longitude, y = Latitude, fill = ..level..,alpha=..level..), bins = 10, geom = "polygon", data = DryCreek) +
scale_fill_gradient(low = "red", high = "green")+
ggtitle("Dry Creek, TN")
Compleated KDE Heatmap
Big Hollow, TN
Big_Hollow <- read_csv("/Users/Kirstenglish/Desktop/Big Hollow.csv")
## Parsed with column specification:
## cols(
## Name = col_character(),
## Latitude = col_double(),
## Longitude = col_double(),
## Number = col_double()
## )
View(Big_Hollow)
Just like before, retrieve a basemap of your site location.
basemap<- get_map(location = "Ridgetop Trail
Ashland City, TN 37015", zoom=15)
## Source : https://maps.googleapis.com/maps/api/staticmap?center=Ridgetop%20Trail%0AAshland%20City,%20TN%2037015&zoom=15&size=640x640&scale=2&maptype=terrain&language=en-EN&key=xxx-txHXwSLjB-UURgX1eLzhNVcsXMRg
## Source : https://maps.googleapis.com/maps/api/geocode/json?address=Ridgetop+Trail%0AAshland+City,+TN+37015&key=xxx-txHXwSLjB-UURgX1eLzhNVcsXMRg
ggmap(basemap)
Put our LOWA encounter cooridinates on the basemap we just created.
#time to put points on the map
BHmap<-ggmap(basemap)+ geom_point(data=Big_Hollow, aes(Longitude, Latitude)) +
geom_point(data = Big_Hollow, aes(Longitude, Latitude), size = 0.1) +
theme(axis.title = element_text(face="bold")) + labs(x="Longitude", y="Latitude") +
theme(panel.border = element_rect(colour = "black", fill=NA, size=1)) + theme_bw()
plot(BHmap)
Put data into spatial points like before to be able to create the KDE and MCP. we will be using CRS (coordinate reference system) which is the format for most GPS data.
Time to make the MCP.
BHmcp <- lapply(BH_sep, FUN = function(x){rgeos::gConvexHull(x)})
##this makes polygon from the list of one when we told R that all variables are from one stream system
BHmcp <- mapply(BHmcp, names(BHmcp),
SIMPLIFY = FALSE,
FUN = function(x,y){x@polygons[[1]]@ID <- y
return(x)})
BHmcp <- do.call(rbind,BHmcp)
BHmcp <- SpatialPolygonsDataFrame(Sr = BHmcp,
data = data.frame(Bird = names(BHmcp)),
match.ID = FALSE)
plot(BHmcp)
Now that we have an MCP for the stream, onto KDE.
bw <- lapply(BH_sep, FUN = function(x){ks::Hlscv(x@coords)})
## Step two: generate kde
BH_kde <-mapply(BH_sep,bw,
SIMPLIFY = FALSE,
FUN = function(x,y){
raster(kde(x@coords,h=y))})
# This code makes a custom function called getContour.
# Inputs:
# kde = kernel density estimate
# prob = probabily - default is 0.95
getContour <- function(kde, prob = 0.95){
# set all values 0 to NA
kde[kde == 0]<-NA
# create a vector of raster values
kde_values <- raster::getValues(kde)
# sort values
sortedValues <- sort(kde_values[!is.na(kde_values)],decreasing = TRUE)
# find cumulative sum up to ith location
sums <- cumsum(as.numeric(sortedValues))
# binary response is value in the probabily zone or not
p <- sum(sums <= prob * sums[length(sums)])
# Set values in raster to 1 or 0
kdeprob <- raster::setValues(kde, kde_values >= sortedValues[p])
# return new kde
return(kdeprob)}
BH_95kde <- lapply(BH_kde,
FUN = getContour,prob = 0.95)
These next plots put MCP on top of KDE
plot(BH_kde[[1]])+
plot(BHmcp[1,],add = TRUE)
## integer(0)
Time to overlap KDE onto the base map. Code found from http://data-analytics.net/cep/Schedule_files/geospatial.html
Compleated KDE Heatmap
Wall Branch, TN
Create the basemap for the site location. This site is Wall Branch stream in Tennesse
basemap<- get_map((location = " 2561 Alex Overlook Way
Clarksville, TN 37043, 36.497818, -87.267428"), zoom=15, maptype = "satellite")
## Source : https://maps.googleapis.com/maps/api/staticmap?center=%202561%20Alex%20Overlook%20Way%0AClarksville,%20TN%2037043,%2036.497818,%20-87.267428&zoom=15&size=640x640&scale=2&maptype=satellite&language=en-EN&key=xxx-txHXwSLjB-UURgX1eLzhNVcsXMRg
## Source : https://maps.googleapis.com/maps/api/geocode/json?address=2561+Alex+Overlook+Way%0AClarksville,+TN+37043,+36.497818,+-87.267428&key=xxx-txHXwSLjB-UURgX1eLzhNVcsXMRg
ggmap(basemap)
LOWA encounter cooridinates on the basemap we just created.
#time to put points on the map
WBmap<-ggmap(basemap)+ geom_point(data=Wall_Branch, aes(Longitude, Latitude)) +
geom_point(data = Wall_Branch, aes(Longitude, Latitude), color= "black", size = 0.1) +
theme(axis.title = element_text(face="bold")) + labs(x="Longitude", y="Latitude") +
theme(panel.border = element_rect(colour = "black", fill=NA, size=1)) + theme_classic()
plot(WBmap)
## Warning: Removed 1 rows containing missing values (geom_point).
## Warning: Removed 1 rows containing missing values (geom_point).
WBpts <- sp::SpatialPoints(coords = cbind(Wall_Branch$Longitude, Wall_Branch$Latitude),
proj4string = sp::CRS("+init=epsg:4326"))
head(WBpts)
## class : SpatialPoints
## features : 1
## extent : -87.26923, -87.26923, 36.49896, 36.49896 (xmin, xmax, ymin, ymax)
## coord. ref. : +init=epsg:4326 +proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0
WB_spdf <- sp::SpatialPointsDataFrame(WBpts, Wall_Branch)
head(WB_spdf)
## # A tibble: 6 x 4
## Name Latitude Longitude Number
## <chr> <dbl> <dbl> <dbl>
## 1 Wall Branch 36.5 -87.3 1
## 2 Wall Branch 36.5 -87.3 2
## 3 Wall Branch 36.5 -87.3 3
## 4 Wall Branch 36.5 -87.3 4
## 5 Wall Branch 36.5 -87.3 5
## 6 Wall Branch 36.5 -87.3 6
WB_sep <- split(x = WB_spdf, f = WB_spdf$Name, drop = FALSE)
Time to make the MCP.
WBmcp <- lapply(WB_sep, FUN = function(x){rgeos::gConvexHull(x)})
##this makes polygon from the list of one when we told R that all variables are from one stream system
WBmcp <- mapply(WBmcp, names(WBmcp),
SIMPLIFY = FALSE,
FUN = function(x,y){x@polygons[[1]]@ID <- y
return(x)})
WBmcp <- do.call(rbind,WBmcp)
WBmcp <- SpatialPolygonsDataFrame(Sr = WBmcp,
data = data.frame(Bird = names(WBmcp)),
match.ID = FALSE)
plot(WBmcp)
Now that we have an MCP for the stream, onto KDE.
bw <- lapply(WB_sep, FUN = function(x){ks::Hlscv(x@coords)})
## Step two: generate kde
WB_kde <-mapply(WB_sep,bw,
SIMPLIFY = FALSE,
FUN = function(x,y){
raster(kde(x@coords,h=y))})
# This code makes a custom function called getContour.
# Inputs:
# kde = kernel density estimate
# prob = probabily - default is 0.95
getContour <- function(kde, prob = 0.95){
# set all values 0 to NA
kde[kde == 0]<-NA
# create a vector of raster values
kde_values <- raster::getValues(kde)
# sort values
sortedValues <- sort(kde_values[!is.na(kde_values)],decreasing = TRUE)
# find cumulative sum up to ith location
sums <- cumsum(as.numeric(sortedValues))
# binary response is value in the probabily zone or not
p <- sum(sums <= prob * sums[length(sums)])
# Set values in raster to 1 or 0
kdeprob <- raster::setValues(kde, kde_values >= sortedValues[p])
# return new kde
return(kdeprob)}
WB_95kde <- lapply(WB_kde,
FUN = getContour,prob = 0.95)
These next plots put MCP on top of KDE
## These next plots put MCP on top of KDE, (territory for this map is very linear)
plot(WB_kde[[1]])+
plot(WBmcp[1,],add = TRUE)
## integer(0)
Time to overlap KDE onto the base map. Code found from http://data-analytics.net/cep/Schedule_files/geospatial.html
WBkde<-WBmap +
stat_density2d(aes(x = Longitude, y = Latitude, fill = ..level..,alpha=..level..), bins = 10, geom = "polygon", data = Wall_Branch) +
scale_fill_gradient(low = "red", high = "green")+
ggtitle("Wall Branch, TN")
Compleated KDE Heatmap
Species Distrubution Map Louisiana Waterthrush
library(dismo)
##
## Attaching package: 'dismo'
## The following object is masked from 'package:ggmap':
##
## geocode
library(rgbif)
library(rdryad)
## Registered S3 method overwritten by 'crul':
## method from
## as.character.form_file httr
library(readxl)
library(ggridges)
##
## Attaching package: 'ggridges'
## The following object is masked from 'package:ggplot2':
##
## scale_discrete_manual
library(viridis)
## Loading required package: viridisLite
library(rasterVis)
## Loading required package: lattice
## Loading required package: latticeExtra
## Loading required package: RColorBrewer
##
## Attaching package: 'latticeExtra'
## The following object is masked from 'package:ggplot2':
##
## layer
extent <- extent(-102,-80,10,45)
LOWA_dismo <- gbif("seiurus", species = "motacilla", ext = extent,
geo = TRUE, sp = TRUE, download = TRUE,
removeZeros = TRUE)
## 938 records found
## 0-300-600-900-938 records downloaded
LOWA_xy <- as.data.frame(cbind(LOWA_dismo@coords[,1],LOWA_dismo@coords[,2]))
colnames(LOWA_xy) <- c("longitude","latitude")
us <- map_data("world")
ggplot(data = LOWA_xy, aes(x=longitude, y=latitude)) +
geom_polygon(data = us, aes(x=long, y = lat, group = group),
fill = "white", color="black") +
geom_point() + xlab("Longitude") + ylab("Latitude") +
coord_fixed(xlim = c(-102,-80), ylim = c(10,45))
bioclim <- getData(name = "worldclim", res = 2.5, var = "bio")
names(bioclim) <- c("Ann Mean Temp","Mean Diurnal Range","Isothermality","Temperature Seasonality","Max Temp Warmest Mo","Min Temp Coldest Mo","Ann Temp Range","Mean Temp Wettest Qtr","Mean Temp Driest Qtr","Mean Temp Warmest Qtr","Mean Temp Coldest Qtr","Annual Precip","Precip Wettest Mo","Precip Driest Mo","Precip Seasonality","Precip Wettest Qtr","Precip Driest Qtr","Precip Warmest Qtr","Precip Coldest Qtr")
bio_extent <- extent(x = c(
min(LOWA_xy$longitude),
max(LOWA_xy$longitude),
min(LOWA_xy$latitude),
max(LOWA_xy$latitude)))
bioclim_extent <- crop(x = bioclim, y = bio_extent)
bioclim_model <- bioclim(x = bioclim_extent, p = LOWA_xy)
presence_model <- dismo::predict(object = bioclim_model,
x = bioclim_extent,
ext = bio_extent)
gplot(presence_model) +
geom_raster(aes(fill=value)) +
geom_polygon(data = us, aes(x= long, y = lat, group = group),
fill = NA, color="black") +
geom_point(data = LOWA_xy, aes(x = longitude, y = latitude), color = "black", alpha = 0.5) +
scale_fill_gradientn(colours=c("brown","yellow","darkgreen"), "Probability") +
coord_fixed(xlim = c(-102,-80), ylim = c(10,45)) +
xlab("Longitude") + ylab("Latitude") + ggtitle("Probability of LOWA Occurrence") +
theme_bw() + theme(plot.title = element_text(hjust = 0.5)) + theme(legend.position = "right") +
theme(panel.grid.major = element_blank(),
panel.grid.minor = element_blank())